home *** CD-ROM | disk | FTP | other *** search
/ Night Owl 6 / Night Owl's Shareware - PDSI-006 - Night Owl Corp (1990).iso / 031a / winlava.zip / SQRT.ASM < prev   
Assembly Source File  |  1991-09-10  |  16KB  |  626 lines
  1. ;---------------------------Module-Header-------------------------------;
  2. ; Module Name: SQRT.ASM
  3. ;
  4. ; Contains FIXED and ULONG square root
  5. ;
  6. ; Created:  Sun 30-Aug-1987 19:28:30
  7. ;
  8. ; NOTE: the 64-bit sqrt code is 386 specific!
  9. ;
  10. ;-----------------------------------------------------------------------;
  11. ;
  12. ;     (C) Copyright Microsoft Corp. 1987-1991.  All rights reserved.
  13. ;
  14. ;     You have a royalty-free right to use, modify, reproduce and 
  15. ;     distribute the Sample Files (and/or any modified version) in 
  16. ;     any way you find useful, provided that you agree that 
  17. ;     Microsoft has no warranty obligations or liability for any 
  18. ;     Sample Application Files which are modified. 
  19. ;
  20.  
  21. ?WIN    = 0
  22. ?PLM    = 0
  23. ?NODATA = 1
  24.  
  25.         .xlist
  26.         include cmacros.inc
  27.         .list
  28.  
  29. UQUAD   struc
  30. uq0     dw      ?
  31. uq1     dw      ?
  32. uq2     dw      ?
  33. uq3     dw      ?
  34. UQUAD    ends
  35.  
  36. PTL    struc
  37. ptl_x    dd    ?
  38. ptl_y    dd    ?
  39. PTL    ends
  40.  
  41. POINTFX struc
  42. x       dd      ?
  43. y       dd      ?
  44. z       dd      ?
  45. w       dd      ?
  46. POINTFX ends
  47.  
  48.  
  49. MAXINT equ 7FFFh
  50. MININT equ 8000h
  51.  
  52.  
  53. ;       The following two equates are just used as shorthand
  54. ;       for the "word ptr" and "byte ptr" overrides.
  55.  
  56. wptr    equ     word ptr
  57. bptr    equ     byte ptr
  58.  
  59. ; The following structure should be used to access high and low
  60. ; words of a DWORD.  This means that "word ptr foo[2]" -> "foo.hi".
  61.  
  62. LONG    struc
  63. lo      dw      ?
  64. hi      dw      ?
  65. LONG    ends
  66.  
  67. FARPOINTER      struc
  68. off     dw      ?
  69. sel     dw      ?
  70. FARPOINTER      ends
  71.  
  72. assert macro   a,b,c,d
  73. endm
  74.  
  75. EAXtoDXAX   macro
  76.         shld    edx,eax,16      ; move HIWORD(eax) to dx
  77.         endm
  78.  
  79. DXAXtoEAX   macro
  80.         ror     eax,16          ; xchg HIWORD(eax) and LOWORD(eax)
  81.         shrd    eax,edx,16      ; move LOWORD(edx) to HIWORD(eax)
  82.         endm
  83.  
  84. .386
  85. sBegin Code
  86.         assumes cs,Code
  87.         assumes es,nothing
  88.         assumes ds,nothing
  89.  
  90. ;---------------------------Private-Macro-------------------------------;
  91. ; set_ov
  92. ;
  93. ; Sets the overflow flag
  94. ;
  95. ; Entry:
  96. ;       reg = word register to use for scratch
  97. ; Returns:
  98. ;       OF set
  99. ; Error Returns:
  100. ;       none
  101. ; Registers Destroyed:
  102. ;       reg
  103. ;-----------------------------------------------------------------------;
  104.  
  105. set_ov  macro   reg
  106.         mov     reg,8000h
  107.         dec     reg
  108.         endm
  109.  
  110. ;---------------------------Public-Routine------------------------------;
  111. ; square_root
  112. ;
  113. ; Computes the Fixed square root of a Fixed.  This algorithm
  114. ; comes from Knuth D.E; Metafont:The Program (Addison-Weseley, 1986)
  115. ; Part 8:Algebraic and Transcendental Functions.
  116. ;
  117. ; Notation used below:
  118. ;
  119. ; Entry:
  120. ;       DX:AX = number to square root
  121. ; Returns:
  122. ;         DX:AX = square root of input
  123. ; Error Returns:
  124. ;       OF = 1
  125. ; Registers Destroyed:
  126. ;       BX,CX
  127. ; Registers Preserved:
  128. ;       DS,ES,SI,DI
  129. ;-----------------------------------------------------------------------;
  130.  
  131.         assumes ds,nothing
  132.         assumes es,nothing
  133.  
  134. cProc   fxsqrt,<PUBLIC,NEAR,NODATA,NONWIN>,<>
  135.         ParmD   fx
  136. cBegin
  137.     mov    ax, fx.lo
  138.         mov     dx, fx.hi
  139.  
  140.         cCall   square_root
  141.  
  142.         mov     dh,dl           ; DX:AX *= sqrt(1000h)
  143.         mov     dl,ah
  144.         mov     ah,al
  145.         xor     al,al
  146. cEnd
  147.  
  148. cProc   ulsqrt,<PUBLIC,NEAR,NODATA,NONWIN>,<>
  149.         ParmD   ul
  150. cBegin
  151.         mov     ax, ul.lo
  152.         mov     dx, ul.hi
  153.  
  154.         cCall   square_root
  155. cEnd
  156.  
  157. ;
  158. ;   compute the square root of EDX:EAX
  159. ;
  160. ;   retult returned in EAX and DX:AX
  161. ;
  162. ;   We dont have a 64-bit square_root function
  163. ;   so approximate it by dividing by 4 until we can use our
  164. ;   32-bit function
  165. ;
  166.  
  167. cProc   sqrt64,<PUBLIC,NEAR,NODATA,NONWIN>,<>
  168. cBegin
  169.         xor     cx,cx
  170.  
  171. test_less_than_32:
  172.         or      edx,edx             ; is it less than 32-bit?
  173.         jz      less_than_32
  174.  
  175.         inc     cx
  176.  
  177.         shrd    eax,edx,2           ; divide by 4
  178.         shr     edx,2
  179.  
  180.         jmp     short test_less_than_32
  181.  
  182. less_than_32:
  183.         ;;
  184.         ;;  EDX:EAX is now less than 32 bits, CX contains count of
  185.         ;;  divisions by 4 we had to do to get here.  call square_root
  186.         ;;  and then re-normalize
  187.         ;;
  188.  
  189.         push    cx
  190.         EAXtoDXAX
  191.         call    square_root
  192.         pop     cx
  193.         jcxz    sqrt64_exit
  194.  
  195. @@:     shl     ax,1
  196.         rcl     dx,1
  197.         loop    @b
  198.  
  199. sqrt64_exit:
  200.         DXAXtoEAX
  201. cEnd
  202.  
  203. cProc   square_root,<PUBLIC,NEAR,NODATA,NONWIN>,<SI,DI>
  204. ;;        ParmD   fx
  205. cBegin
  206. ;;        mov     ax, fx.lo
  207. ;;        mov     dx, fx.hi
  208. ; check for < 1 since this algorithm returns 0000:0001
  209.  
  210.         or      dx,dx
  211.         jnz     non_zero_arg
  212.         cmp     ax,1
  213.         ja      non_zero_arg
  214.         jmp     square_root_exit
  215.  
  216. negative_arg:                           ; can't have number negative
  217.         xor     ax,ax
  218.         cwd
  219.         set_ov  bx
  220.         jmp     square_root_exit
  221.  
  222. non_zero_arg:
  223.         cCall   ulNormalize
  224. ;;;;;;;;jcxz    negative_arg
  225.         jz      exit_relay
  226.         shr     cx,1
  227.         jc      add_shifts_only
  228.         shr     dx,1
  229.         rcr     ax,1
  230.         dec     cx
  231. add_shifts_only:
  232.  
  233. ;;
  234. ;; use 23 for a FIXED square_root, and 15 for a ULONG square_root!
  235. ;;
  236.  
  237. if 0
  238.         mov     ch,23                   ; FIXED square_root
  239. else
  240.         mov     ch,15                   ; ULONG square_root
  241. endif
  242.  
  243.         sub     ch,cl
  244.         mov     cl,8
  245.         sub     ch,cl
  246.         jg      more_than_8             ; CL = Min(count,8)
  247.         add     cl,ch
  248.         xor     ch,ch                   ; CH = Max(count-8,0)
  249. more_than_8:
  250.  
  251.  
  252.         xchg    ax,si                   ; DI:SI = x
  253.         mov     di,dx
  254.  
  255.         xor     ax,ax                   ; AX = y = lower bound = 0
  256.  
  257.         mov     bx,2                    ; BX = q = estimate of root = 2
  258.  
  259.         add     si,si                   ; intiialize y = top bit of x
  260.         adc     di,di
  261.         adc     ax,ax
  262.  
  263. first_8_loop:
  264.         add     si,si                   ; move 2 bits from top
  265.         adc     di,di                   ; of x to bottom of y
  266.         adc     ax,ax
  267.         assert  NC
  268.  
  269.         add     si,si
  270.         adc     di,di
  271.         adc     ax,ax
  272.         assert  NC
  273.  
  274.         sub     ax,bx                   ; AX = y - q
  275.         jle     y_neg_or_zero_8
  276.  
  277. ;!!!CR loops can be cleaned up some
  278.  
  279.         add     bx,bx                   ; BX = 2q
  280.         sub     ax,bx                   ; AX = y - 3q
  281.         jg      y_greater_q_8
  282.  
  283.         add     ax,bx                   ; AX = y - q
  284.  
  285.         dec     cl
  286.         jnz     first_8_loop
  287.  
  288.         jcxz    first_8_was_enough
  289.         jmp     short done_with_1st_8
  290.  
  291. y_greater_q_8:
  292.         add     bx,2                    ; BX = 2q + 2 (or 2q if jg not taken)
  293.  
  294.         dec     cl
  295.         jnz     first_8_loop
  296.  
  297.         or      ch,ch
  298.         jnz     done_with_1st_8
  299.  
  300. first_8_was_enough:
  301.         xchg    ax,bx                   ; DX:AX = q
  302.         xor     dx,dx
  303.         shr     ax,1                    ; root = q >> 1
  304.         or      ax,ax                   ; clear overflow flag
  305. exit_relay:
  306.         jmp     square_root_exit
  307.  
  308. y_neg_or_zero_8:                       ; come here if y <= 0
  309.         dec     bx
  310.         add     bx,bx                   ; BX = 2q - 2
  311.         add     ax,bx                   ; AX = y + q - 2
  312.  
  313.         dec     cl
  314.         jnz     first_8_loop
  315.         jcxz    first_8_was_enough
  316.  
  317. ;si is now empty, di contains all signifigant bits
  318. done_with_1st_8:
  319.  
  320.         mov     cl,4
  321.         sub     ch,cl
  322.         jg      second_4_loop           ; CL = Min(count,4)
  323.         add     cl,ch
  324.         xor     ch,ch                   ; CH = Max(count-8,0)
  325.  
  326. second_4_loop:
  327.                                         ; move 2 bits from top
  328.         add     di,di                   ; of x to bottom of y
  329.         adc     ax,ax
  330.         assert  NC
  331.  
  332.         add     di,di
  333.         adc     ax,ax
  334.         assert  NC
  335.  
  336.         sub     ax,bx                   ; AX = y - q
  337.         jle     y_neg_or_zero_2nd_4
  338.  
  339.         add     bx,bx                   ; BX = 2q
  340.  
  341. ;!!!CR use compare and jump rather than sub, branch, and restore
  342.         sub     ax,bx                   ; AX = y - 3q
  343.         jg      y_greater_q_2nd_4
  344.  
  345.         add     ax,bx                   ; AX = y - q
  346.  
  347.         dec     cl
  348.         jnz     second_4_loop
  349.  
  350.         jcxz    second_4_was_enough
  351.         jmp     short done_with_2nd_4
  352.  
  353. y_greater_q_2nd_4:
  354.         add     bx,2                    ; BX = 2q + 2 (or 2q if jg not taken)
  355.  
  356.         dec     cl
  357.         jnz     second_4_loop
  358.  
  359.         or      ch,ch
  360.         jnz     done_with_2nd_4
  361.  
  362.  
  363. second_4_was_enough:
  364.         xchg    ax,bx                   ; DX:AX = q
  365.         xor     dx,dx
  366.         shr     ax,1                    ; root = q >> 1
  367.         or      ax,ax                   ; clear overflow flag
  368.         jmp     square_root_exit
  369.  
  370.  
  371. y_neg_or_zero_2nd_4:                       ; come here if y <= 0
  372.         dec     bx
  373.         add     bx,bx                   ; BX = 2q - 2
  374.         add     ax,bx                   ; AX = y + q - 2
  375.  
  376.         dec     cl
  377.         jnz     second_4_loop
  378.         jcxz    second_4_was_enough
  379.  
  380. done_with_2nd_4:
  381.  
  382.         xor     dx,dx                   ; DX:AX = y   SI:BX = q
  383.  
  384.         mov     cl,4
  385.         sub     ch,cl
  386.         jg      third_4_loop             ; CL = Min(count,4)
  387.         add     cl,ch
  388.         xor     ch,ch                   ; CH = Max(count-13,0)
  389.  
  390. third_4_loop:
  391. ;!!!CR Add comment noting that 32 bit math is being used now
  392.                                         ; move 2 bits from top
  393.         add     di,di                   ; of x to bottom of y
  394.         adc     ax,ax
  395.         adc     dx,dx
  396.         assert  NC
  397.  
  398.         add     di,di
  399.         adc     ax,ax
  400.         adc     dx,dx
  401.         assert  NC
  402.  
  403.         sub     ax,bx                   ; DX:AX = y - q
  404.         sbb     dx,si
  405.         js      y_negative_3rd_4
  406.         jz      y_might_be_zero_3rd_4
  407.  
  408. sorry_y_not_zero_3rd_4:
  409.         add     bx,bx
  410.         adc     si,si                   ; SI:BX = 2q
  411.         sub     ax,bx                   ; DX:AX = y - 3q
  412.         sbb     dx,si
  413.         jg      y_greater_q_3rd_4
  414.         jl      y_less_or_equal_q_3rd_4
  415.         or      ax,ax
  416.         jnz     y_greater_q_3rd_4
  417.  
  418. y_less_or_equal_q_3rd_4:
  419.         add     ax,bx                   ; DX:AX = y - q
  420.         adc     dx,si
  421.  
  422.         dec     cl
  423.         jnz     third_4_loop
  424.  
  425.         jcxz    third_4_was_enough
  426.         jmp     short done_with_3rd_4
  427.  
  428. y_greater_q_3rd_4:
  429.         add     bx,2
  430.         adc     si,0
  431.  
  432.         dec     cl
  433.         jnz     third_4_loop
  434.  
  435.         or      ch,ch
  436.         jnz     done_with_3rd_4
  437.  
  438. third_4_was_enough:
  439.         xchg    ax,bx                   ; DX:AX = q
  440.         mov     dx,si
  441.         shr     dx,1                    ; root = q >> 1
  442.         rcr     ax,1
  443.         or      ax,ax                   ; clear overflow flag
  444.     jmp    square_root_exit ;;; short!
  445.  
  446. y_might_be_zero_3rd_4:
  447.         or      ax,ax
  448.         jnz     sorry_y_not_zero_3rd_4
  449.  
  450. y_negative_3rd_4:                       ; come here if y <= 0
  451.  
  452.         sub     bx,1
  453.         sbb     si,0
  454.         add     bx,bx                   ; DL:BX = 2q - 2
  455.         adc     si,si
  456.  
  457.         add     ax,bx
  458.         adc     dx,si                   ; DH:AX = y + q - 2
  459.  
  460.         dec     cl
  461.         jnz     third_4_loop
  462.         jcxz    third_4_was_enough
  463.  
  464.  
  465. done_with_3rd_4:
  466.  
  467.         xchg    ch,cl
  468.  
  469. last_7_loop:
  470.         add     ax,ax
  471.         adc     dx,dx
  472.  
  473.         add     ax,ax
  474.         adc     dx,dx
  475.  
  476.         sub     ax,bx                   ; DX:AX = y - q
  477.         sbb     dx,si
  478.         js      y_negative_last_7
  479.         jz      y_might_be_zero_last_7
  480.  
  481. ;!!!CR Clean up the exit so that a single exit point is used.
  482.  
  483. sorry_y_not_zero_last_7:
  484.         add     bx,bx
  485.         adc     si,si                   ; SI:BX = 2q
  486.         sub     ax,bx                   ; DX:AX = y - 3q
  487.         sbb     dx,si
  488.         jg      y_greater_q_last_7
  489.         jl      y_less_or_equal_q_last_7
  490.         or      ax,ax
  491.         jnz     y_greater_q_last_7
  492.  
  493. y_less_or_equal_q_last_7:
  494.         add     ax,bx                   ; DX:AX = y - q
  495.         adc     dx,si
  496.  
  497.         loop    last_7_loop
  498.  
  499.         xchg    ax,bx                   ; DX:AX = q
  500.         mov     dx,si
  501.         shr     dx,1                    ; root = q >> 1
  502.         rcr     ax,1
  503.         or      ax,ax                   ; clear overflow flag
  504.         jmp     short square_root_exit
  505.  
  506. y_greater_q_last_7:
  507.         add     bx,2
  508.         adc     si,0
  509.  
  510.         loop    last_7_loop
  511.  
  512.         xchg    ax,bx                   ; DX:AX = q
  513.         mov     dx,si
  514.         shr     dx,1                    ; root = q >> 1
  515.         rcr     ax,1
  516.         or      ax,ax                   ; clear overflow flag
  517.         jmp     short square_root_exit
  518.  
  519. y_might_be_zero_last_7:
  520.         or      ax,ax
  521.         jnz     sorry_y_not_zero_last_7
  522.  
  523. y_negative_last_7:                       ; come here if y <= 0
  524.  
  525.         sub     bx,1
  526.         sbb     si,0
  527.         add     bx,bx                   ; DL:BX = 2q - 2
  528.         adc     si,si
  529.  
  530.         add     ax,bx
  531.         adc     dx,si                   ; DH:AX = y + q - 2
  532.  
  533.         loop    last_7_loop
  534.  
  535.         xchg    ax,bx                   ; DX:AX = q
  536.         mov     dx,si
  537.         shr     dx,1                    ; root = q >> 1
  538.         rcr     ax,1
  539.  
  540.         or      ax,ax                   ; clear overflow flag
  541. square_root_exit:
  542. cEnd
  543.  
  544. ;---------------------------Public-Routine------------------------------;
  545. ; ulNormalize
  546. ;
  547. ; Normalizes a ULONG so that the highest order bit is 1.  Returns the
  548. ; number of shifts done.  Also returns ZF=1 if the ULONG was zero.
  549. ;
  550. ; Entry:
  551. ;       DX:AX = ULONG
  552. ; Returns:
  553. ;       DX:AX = normalized ULONG
  554. ;       CX    = shift count
  555. ;       ZF    = 1 if the ULONG is zero, 0 otherwise
  556. ; Registers Destroyed:
  557. ;       none
  558. ;-----------------------------------------------------------------------;
  559.  
  560.         assumes ds,nothing
  561.         assumes es,nothing
  562.  
  563. cProc   ulNormalize,<PUBLIC,NEAR>
  564. cBegin
  565.  
  566. ; shift by words
  567.  
  568.         xor     cx,cx
  569.         or      dx,dx
  570.         js      ulNormalize_exit
  571.         jnz     top_word_ok
  572.         xchg    ax,dx
  573.         or      dx,dx
  574.         jz      ulNormalize_exit        ; the zero exit
  575.         mov     cl,16
  576.         js      ulNormalize_exit
  577. top_word_ok:
  578.  
  579. ; shift by bytes
  580.  
  581.         or      dh,dh
  582.         jnz     top_byte_ok
  583.         xchg    dh,dl
  584.         xchg    dl,ah
  585.         xchg    ah,al
  586.         add     cl,8
  587.         or      dh,dh
  588.         js      ulNormalize_exit
  589. top_byte_ok:
  590.  
  591. ; do the rest by bits
  592.  
  593.         inc     cx
  594.         add     ax,ax
  595.         adc     dx,dx
  596.         js      ulNormalize_exit
  597.         inc     cx
  598.         add     ax,ax
  599.         adc     dx,dx
  600.         js      ulNormalize_exit
  601.         inc     cx
  602.         add     ax,ax
  603.         adc     dx,dx
  604.         js      ulNormalize_exit
  605.         inc     cx
  606.         add     ax,ax
  607.         adc     dx,dx
  608.         js      ulNormalize_exit
  609.         inc     cx
  610.         add     ax,ax
  611.         adc     dx,dx
  612.         js      ulNormalize_exit
  613.         inc     cx
  614.         add     ax,ax
  615.         adc     dx,dx
  616.         js      ulNormalize_exit
  617.         inc     cx
  618.         add     ax,ax
  619.         adc     dx,dx
  620. ulNormalize_exit:
  621. cEnd
  622.  
  623. sEnd    Code
  624.  
  625.         end
  626.